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ABSTRACT 



As pan of the international Greenland Sea Project, Woods Hole Oceanographic 
Institute and Scripps Institution of Oceanography deployed a six transceiver ocean 
acoustic tomography array to monitor ocean ventilation and circulation over the 
I9SS-89 winter cooling season. A stochastic inverse method computer code which at- 
tains a solution by mininrizing mean square error is used to perform inversions of the 
Greenland Sea tomography data. 

A computer simulated ocean is used to evaluate various aspects of system perform- 
ance. W'c first consider the advantages and problems associated with using a ray theory 
based algorithm to establish the fonvard problem for the Greenland Sea tomography 
array. Next, we made two adjustments to our inversion code and discuss the eflects on 
system performance. The first adjustment allows for layers of dilTerent thicknesses in the 
inverse solution to increase the density of estimates in regions of interest. The second 
adjustment allows the estimator to expect variability of the unknown field to decrease 
exponentially with depth. 

Our results show the ray theory based algorithm is an adequate method of modeling 
ray paths in the Greenland Sea, but has limitations. Reliability of ray paths degrades 
as launch angles become s.iallower and if strong gradients and rapidly changing gradi- 
ents in sound speed are present in the vicinity of the transceiver elements. We also find 
varxing the thickness of layers in the solution allows us to examine the more variable 
upper ocean in greater detail without increasing computational efibrt. However, this 
adjustment alone has the undesirable side effect of shifting system resolution towards the 
lower ocean. This shift in resolution is oflset by informing the estimator about the ver- 
tical variability distribution of the unknown field. This a priori knowledge is 
parameterized by the covariance function of the unknown field. Uncertainty in knowing 
the true variability distribution aflects model performance. The inverse solution is more 
sensitive to underestimating than overestimating the true value of folding depth. The 
model is also more sensitive to both underestimation and overestimation at small true 
folding depths. 

A set of Greenland Sea data between one tranceiver pair was processed by Woods 
Hole Oceanographic Institute. Although only three groups of eigenrays are involved, 
initial inversion results indicate the estimator detects seasonal changes and synoptic scale 
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events occurring at lime scales greater than 20 days, however, solutions show wide fluc- 
tuations at time scales shorter than 20 days. 
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I. INTRODUCTION 



A. OCEAN ACOUSTIC TOMOGRAPHY 

Tomography is a method to remotely sense interior structures. It uses 
electromagnetic, seismic, or sound waves to propagate through and probe relatively 
transparent media. For e.xample, x-rays are commonly used in medical computer as- 
sisted tomography (C.AT) scans and seisntic waves induced by surface sources are used 
to study the interior of the earth (Backus and Gilbert, 1967). Using analagous methods, 
Munk and Wunsch (1979) introduced ocean acoustic tomography as a means of moni- 
toring mesoscale fluctuations in ocean basins. Sound waves, which arc sensitive to 
density changes and currents but transit well through the ocean, gather information as 
they encounter structures in the ocean interior. The data is received in the form of a 
peturbated travel time from a source to a receiver. Using inverse techniques, the best 
estimate as to what structure could result in the observed data is constructed. 

Ocean acoustic tomography has several practical advantages over conventional 
hydrographic study techniques, as cited by Chiu ei al. (19S7). The monitoring system 
can be established as a sentipermanent, continuous observing system which is not greatly 
alFected by weather. It has high temporal resolution and can cover an extensive volume 
of ocean simultaneously with only a few moorings, which reduces deployment and 
maintenance costs. Another advantage stated by Munk and Wunsch (1979) is that tlie 
amount of additional information gained by each additional mooring increases more 
rapidly than that gained by additional "spot" measurements since each mooring sets up 
distinct new paths to each of the previously set moorings. The gain in spot measure- 
ments is more or less linear with the number of instruments deployed. 

Ocean acoustic tomography has been successfully demonstrated in a variety of ap- 
plications since its introduction. Some examples include mapping mesocale midocean 
eddy fields (Cornuelle ei al.. 19S5), analysis of planetarx' waves (Chiu ei «/.,19S7), 
studying of ocean currents (DeFerrari ei al., 19S6), and surface wave analysis (Uynch ci 
n/..19S7). It has also been chosen as a monitoring tool as part of the Greenland Sea 
Project. Its application in the Greenland Sea is the impetus for this thesis. 

The Greenland Sea is a well known contributor to world ocean ventilation. Several 
conceptual theories have been developed to describe the process, but direct observations 
of ventilation are scarce. The Greenland Sea Project was put forward by the Arctic 
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Ocean Sciences Board to gain understanding of the processes relevant to deep con- 
vection (Meincke.l9S9). 

B. THE GREENLAND SEA PROJECT 

The Greenland Sea Project (GSP) is a five year (1987-1992), international scientific 
study to understand the large scale, long term interactions of air, sea, and ice in the 
Greenland Sea. 1 he five key elements of the GSP are: (1) a study of air-sea-ice inter- 
action. (2) a study of ocean ventilation. (3) a study of ocean circulation and mixing, (4) 
a study of atmospheric energetics, and (5) a study of biological processes. The role of 
ocean acoustic tomography is to assist in the studies of ocean ventilation and circu- 
lation. The array of transceiver moorings is designed to measure changes in the inte- 
grated properties of the Greenland Sea gyre through a winter cooling season and may 
also provide valuable information in measuring the response of the wind driven gyre to 
changes in the curl of the wind stress (Greenland Sea Science Planning Group. 1986). 

,'\ six transceiver array was successfully deployed in the Greenland Sea during 
September-October, 1988. Figure 1 depicts the planned deployment area and geometiw 
of the array. The moorings were shifted slightly from their original planned location to 
accommodate for the rough bathymetry found at the original .Mooring 3 site. Mooring 
2 was redeployed due to leakage of some O-ring seals. The redeployed mooring is des- 
ignated here as .Vlooring 2a. '1 he geodic coordinates of the moorings, the source depths, 
and the receiver depths are given in Table 1. (Worchester and Howe. 1989). 

C. THESIS OBJECTIVES 

There are three basic objectives of this thesis. The first is to discuss the advantages 
and problems of predicting ray paths between array transceivers using an algoritm based 
on ray theory. I'he Greenland Sea oflers an acoustic environment which, in conjunction 
with the array geometry, makes ray path determination challenging to model. Accurate 
ray path prediction is an essential element of establishing the "fonvard" problem. 
Chapter III contains our discussion on this topic. 

The second objective is to develop an inversion code appropriate for application to 
Greenland Sea tomography data. We begin with an inversion code originated by Chiu 
and modified by Kao (1989). The code, or estimator, treats the ocean as a volume which 
has been subdivided into discrete boxes. The boxes have equal horizontal dimensions 
and the vertical layers of the ocean are equally spaced. An estimate of the unknown 
variable is calculated for each box. The estimator also treats the ocean as a statistically 




Figure 1. Planned Dcployincnt Site of llie Greenland Sea Tomography Array 
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Table !. TRANSCEI\’ER LOCATIONS 



-Mooring 


Lat-Long Position 


X-Y Position (km) 


Source 
Depth (m’ 


Receiver 
Depth (m 


1 


75'58.0S'.V.OOr5O.00Tr 


154.790 220.317 


99.7 


150.4 


2 


75°03.69'.V.O00M0.25'£ 


224.683 119.080 


94.0 


144.7 


2a 


75'03.88'.V.OOO'38.20’£ 


223.715 119.491 


97.4 


148.1 


3 


74'09.38'.V.00r52.90Tr 


148.972 18.146 


96.0 


146.7 


4 


74'28.90'.V.O05'47.30Tr 


33.270 60.958 


94.6 


117.7 


5 


75'34.27'A’.O06°O7.70Tr 


34.547 182.902 


101.5 


124.6 


6 


75'03.6O'.V.0O2°5S.00Tr 


1 2O.000 1 20.000 


92.5 


143.2 



homogeneous volume. Therefore it expects variability (RMS values of the unknown 
variable) in any region of the ocean to be the same. 

We have made two fundamental adjustments to the inversion code and discuss how 
they eflcct system performance. Our first adjustment is to vary the layer thickness so 
that there is a greater density of boxes in the region of interest. This adjustment allows 
for more detail in the region of interest without increasing the computational effort. 
Our second adjustment tells the estimator to expect a depth dependent variability dis- 
tribution of the unknown variable. In this way. the estimator produces solutions which 
arc more realistic statistically. The efl'ects of the adjustments on system performance is 
evaluated by conducting inversions on computer simulated oceans with known statistical 
parameters. The results of these studies are also contained in Chapter III. 

For convenience in coding, we use a Cartesian coordinate system. Based on ranges 
computed for the WGS84 spheroid (Worchester and Howe, 1989). we converted the 
latitude and longitude positions to an X'i' plane which closely approximates the geodic 
orientation. The plane is 240 km by 240 km and centered on Mooring 6. The XY co- 
ordinates are included in Table 1. Though there is some slight distortion in the overall 
shape of the array, the error between the geodic ranges and the model ranges is less than 
15 meters along each path, as shown in Table 2 below. 
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Table 2. MOORING RANGES (KM) AND RANGE ERRORS (M): Ranges 

(above diagonal) are based on the WGSS4 spheroid, range errors (below 
diagonal) are ranges minus spheroid ranges. 



■M coring 


1 


2 


2a 


3 


4 


5 


6 


1 


0 


123.033 


122.120 


202.258 


200.405 


125.922 


106.179 


2 


-13 


0 


1 .039 


126.177 


200.057 


200.576 


104.701 


2a 


13 


13 


0 


125.927 


199.232 


199.500 


103.706 


3 


-3 


-3 


-1 


0 


123.367 


200.594 


105.896 


4 


1 


-14 


5 


2 


0 


121.947 


104.921 


5 


8 


-14 


13 


-1 


4 


0 


106.109 


6 


-1 


-14 


10 




_2 


-1 


0 



Our third objective is to test the inversion code with observed data derived from the 
Greenland Sea tomographic array. The studies on system performance discussed earlier 
are conducted using computer simulated oceans with known characteristics. However, 
the simulated oceans arc limited in the scope to which they can mimmick the actual 
conditions found in the Greenland Sea. Using actual Greenland Sea travel time data is 
the only way to determine if the model response will be realistic. Currently there is a set 
of data derived from signals transmitted between Mooring 4 and Mooring 5 available for 
input into our model. The data set was processed by Woods Hole Oceanographic In- 
stitute (^^'HO^). Though this will not enable us to form significant conclusions about 
the array site from such a small data base, ingesting real data into our model is helpful 
in preparing the inversion code for the complete data set when it becomes available. 
Our results arc discussed in Chapter IV. 

To complete our study, we provide the following information. In Chapter II, we 
present an overview of the basis on which our system is developed. Here, the math- 
ematical foundation of the tomographic "forward" and "inverse" problems is discussed. 
Then in Chapter V. we summarize our final conclusions and make recommendations for 
future improvements related to this work. 
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II. OVERVIEW OF ACOUSTIC TOMOGRAPHY PRINCIPLES 



In general, the application of tomogrgaphy to monitoring the ocean interior is done 
in two parts. In the "forward" problem, we establish a model describing the physical 
relationship between the observed data and the unknown variables which are to be esti- 
mated. We will cast the forward problem as a Fredholm Integral Equation of the First 
Kind in the form 



di=\ g{s)J[s)ds ->r ei /=l,2,...,m (2.1) 

''path 

where d, is the data accumulated along the path. g,(s) is the data kernel which contains 
the physics and operates on the unknown, J{s) is the unknown structure which is to be 
estimated, and e, is the experimental noise accumulated along the /''■ path. 

The "inverse" problem is to deternrine a solution to the unknown structure which fits 
the data. An undetermined problem, one generally finds many solutions that fit the data 
and must choose an optimal solution based on some objective criteria. For the 
Greenland Sea transceiver array, the observed data arc travel times of sound energy 
along reciprocal acoustic multipaths which "connect" sources and receivers. From these 
data, we can estimate both sound speed perturbations with respect to a reference ocean 
(density tomography), and ocean currents (current tomography). 

A. THE FORWARD PROBLEM 

Sound propagation in the ocean is a function of temperature, pressure, and salinity 
as well as a function of the current field. For the forward problem, we need to determine 
the acoustic multipaths which sound energy follows from source to receiver. To model 
these "eigenrays", we use an algorithm based on ray theory. Ray theory provides a 
simple physical description of the acoustic multipaths and the modeling equations are 
straightforward. We start with the eikonal equation, 

1 V<iy{x^-,z) 1^ = «(x,v.r)‘ (2.2) 

where O is the acoustic phase defining a wave front and n is the index of refraction. 
Following the example of Ziomek (1985), we neglect horizontal refraction so that the 
predicted ray paths lie on a vertical plane passing through the source and receiver. 
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Following the derivation used by Kao (1989), and treating the reference ocean as a 
motionless medium, we arrive at the key equation for modeling ray paths in the 
Greenland Sea. 



(*R 






' »( 2 ) 
COS 



s 2 



- 1 dR 



(2.3) 



where z is the depth of the ray at range R, 6 is the ray path angle with respect to the 
horizontal, and the subscript o indicates conditions at the acoustic source. 

The speed of the wave front through the true ocean is affected by both the sound 
speed and current structure along the path. We express this as 

= C(.T,V.2./) + V(.YJ-.2,/) . s (2.-4) 

where c is sound speed. V is the ocean current field, and i is a unit vector tangent to the 
ray path, s . and pointed in the direction of sound energy propagation. From Equation 
(2.4). we can determine the sound travel time along an eigenray path. 




We further separate c into the sum of sound speed in a reference ocean and sound speed 
perturbations with respect to this reference, 



c = c^ + Sc (2.6) 

Travel time can also be separated into the sum of travel time along an eigenray in 
a reference ocean and the perturbation in travel lime due to effects in the ocean not re- 
presented by the reference conditions. 



T= T^ + ST 



(2.7) 



Substituting Equations (2.6) and (2.7) into Equation (2.5), we have 
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7 . + *^^ = 



path 






( 2 . 8 ) 



Under the mild restriction that | dc + V . i | < < | cj .we can neglect the higher order 
terms of a Taylor series expansion of the term in parenthesis and arrive at 



In density tomography, which uses one way travel times, the effect of ocean currents 
is much smaller than the elfects of temperature, pressure, and salinity, i.e. 



In current tomography, we must consider reciprocal paths. The effects of temper- 
ature. pressure, and salinity are the same in either direction whereas the efl'ects of the 
current field will be opposite under these circumstances. Therefore, imbedded informa- 
tion about the currents may be extracted by subtracting the two reciprocal path 
equations. 

Consider the the travel time for sound along a path indicated by =s> calculated from 
Hquation (2.9), 




path 



(2.9) 



I • i I < < I (5c 1 . Since we have defined = 
linearized form 



Equation (2.9) reduces to the 




( 2 . 10 ) 




( 2 . 11 ) 



and for its reciprocal path indicated by 



S 



( 2 . 12 ) 



To + ^l = 



^paih 



^ I- 



Sc + V • s' 



(is 



Since V . i“ = — V • j-* at all points on the paths, subtracting Equation (2.12) from 
Equation (2.1 1 ) yields a linearized expression for current, 



•''paili 

Thus far. the forward model has treated travel time data only as noise free meas- 
urements. In practice, experimental noise corrupts the signal and contributes to the er- 
ror in the estimates. Factors such as tidal forces which cause displacement of transceiver 
elements and internal waves are sources of experimental noise. Signal processing is used 
to minimize these efTects. Details of signal processing techniques used for ocean 
tomography can be found in Spindel (1985). 

Including experimental noise in our model, we now have expressions which relate 
the unknown variables to practical, observable data. For density tomography, we have 



CV . 



(2.13) 



6T = 




)((h:)<fs- -f 



and for current tomography, we have 



(2.1d) 



(V.s)i/s + e- (2.15) 

The similarity of the data kernels, which contain the physics of the model, shows the 
close relationship between the two areas of tomography. 
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B. THE INT ERSE PROBLEM 



Once the forward problem has been established, we search for unknown structures 
which fit the observed data. Since we cannot sufTiciently sample the ocean to have a 
complete set of independent data, the problem is underdetermined, with an infinite 
number of possible solutions that will fit the data. To reach a unique solution, we must 
devise an estimator which imposes additional constraints on the system. Choosing an 
estimator based on the Gauss-.MarkolT Theorem allows the use of statistical information 
to evaluate system performance. 

1. The Estimator 

We approximate the continuous system in a discretized form. This allows us to 
cast the problem into matrix algebra for straightforward implementation of the estimator 
in a numerical computer. Applying the Gauss-Markoff Theorem (Liebelt,1967j and 
following the works of Cornuelle ci al. (1985) and Chiu ei al. (1987), our criterion for the 
"best" solution of the unknown structure is one that is linear with the data and has 
minimal mean square error with respect to the true solution. Additional contraints on 
the system are imposed through the specification of an a priori covariance matrix. 

Posing the forward problem in discretized form, we have 

d^\f+e (2.16) 

where the m dimensional column vector, t/. contains the data derived from travel times 
of m resolvable eigenrays. The /? dimensional column vector,/, is the unknown structure 
to be estimated. The ocean has been divided into n discrete boxes, within each the un- 
known variable. Sc or V. is assumed to be constant. In our treatment, the discrete boxes 
have equal dimensions in the horizontal plane, forming squares. However, the vertical 
dimension varies with depth, forming uniform layers of different thicknesses. The ex- 
perimental noise, e, which contaminates the data along each ray path, is also m dimen- 
sional. The linear operator matrix, A. has m x « dimensions. The rows of the operator 
matrix are analagous to a set of data kernels of the continuous case. They contain the 
physics which relates the data and unknown fields. 

From the Gauss-Markoff Theorem, the optimal estimate of an unknown pa- 
rameter can be calculated from 



/ = c,A^c;‘j/ 

where the covariance matrix of the total error, e. of the estimate is 



(2.17) 
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( 2 . 18 ) 



C, = - (C^a"')(AC^a''+ C,)-’(CyA"')^ 



The diagonal elements of C, are the minimum mean square errors of the estimate. C, is 
the covariance matrix of experimental noise, which we assume is known, and Cf is the 
covariance matrix of the unknown variable. It is through that we insert a priori in- 
formation in the estimator. To reach an optimal solution, we assume the unknown fields 
have Gaussian distributions and are statistically homogeneous throughout each layer. 
We express these conditions through a covariance function, Q , which is cast into matrix 
form as . The covariance function is expressed as. 



C^= o{x',y', . 



:')o(x",y", z") expj 




where the separation of two points in each dimension is given as 
Ax = .v" — .v'. Ay=y" — y', and Az = z" — z' , respectively. The dimensional correlational 
lengths. L,. and L, tell the estimator how well information at one point relates to 
surrounding points in the system. It is optimal to use the true correlation lengths, but 
the true values may not be well known. To avoid overconstraining the system, i.e. to 
avoid correlating information which is not truly correlated, we can use conservatively 
small values of L,. L^. and L,. We characterize the vertical variability distribution of the 
unknown field through a, which is treated by the estimator as an exponentially decaying 
function with depth, parameterized by a surface value and an e -folding depth. In each 
layer, an average value based on this curve is used. With a priori knowledge in place, 
the inversion code can generate a unique least mean square error solution from observed 
data as well as error bars and system resolution measures. 

2. Measurement of Error 

The total error in the estimate can be cast as two statistically independent terms; 
bias, b = </> —f. and random error, Af =f — <f > . The covariance of error is simply 
the sum of its components, 

C, = <^^> + C^; (2.20) 



where 

C^; = C,A^C;’ AC, • (2.21) 
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Following Chiu et al. (19S7) and Wiggins (1972). we apply a powerful tool called singular 
value decoinposition of matrices. The linear operator matrix. A, is scaled and factorized 
into a diagonal matrix of eigenvalues, A, and two matrices of associated eigenvectors. 
U and V such that. 

_ j_ j_ 

C, 2 xc/ = U A (2.22) 

Each diagonal element of A. has two associated eigenvectors, u, and v, , which are the 
C columns of U and V, respectively. The eigenvectors within U and V are sets of 
orthogonal basis vectors in the data space and unknown parameter space, respectively. 

.Applying matrix decomposition as in Chiu ei al. (I9S7), we now write the 
covariances of total error and random error as 

X ± 

= c/[l - VA(I + A^)"‘AVT C/ (2.23) 

and 

X _L 

Ca/ = + A"r'A^(l + aV'vT C/ (2.24) 

where I is an « x w identity matrix. We derive an expression for bias by taking the ex- 
treme case where all eigenvalues approach infinity so that C^; aproaches zero. From 
Equation (2.23) we have, 

X ± 

< bb'^ > = Cf\l - (2.25) 



3. Measurement of Resolution 

Besides a measure of error, we also evaluate system performance through its 
ability to resolve ocean features. Continuing wdth singular value decomposition of ma- 
trices, we obtain an expression for the optimal estimate in terms of eigenvalues, 



C/ 





(2.26) 



where k is the number of non-zero eigenvalues. Equation (2.26) shows the solution is a 
weighted sum of eigenvectors where the eigenvalues control the weighting. The ex- 
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pression in parenthesis represents the components of the data vector which has been 
expanded into a linear combination of weighted eigenvectors. As demonstrated by Chiu 
et al. (19S7). the /? values are analagous to signal-to-noise ratios for each component 
of tlie expanded data vector. For /f > > I, signal dominates noise in that particular 
component, and the reverse is true for < < 1. From Equation (2.26) we see that the 
effect of noise on the estimate is ntinimized by downweighting those components asso- 
ciated with high noise levels. 

We now define the resolution matrix, R, which has w x « dimensions, 

R = VA(1 -t- aV'AV^ (2.27) 

The /•■'' row of R is the resolution kernel of the box in the discretized ocean. Substi- 
tuting into Equation (2.23). we get 

_L j_ 

C, = - R)C^2 (2.28) 



Ideally, if the resolution matrix is an identity matrix, the system has no error. In the 
limit, as we approach the continuous case with complete and noise free data, the resol- 
ution kernels approach delta functions and the system would be perfectly resolved. In 
practice, however, the resolution kernels have side lobes and amplitudes less than unity 
along the diagonal. In Equation (2. 28), the relationship between resolution and error is 
linear. 1 herefore, we may use the resolution kernel peak value as a simple measure of 
local resolution. 

.Another measure of system performance is a measure of the size of features that 
can be adequately resolved by the system. Chiu ei al. (1987) has defined scale meas- 
urements called the horizontal and vertical minimum resolution lengths, H and V, re- 
spectively. Each represents the distance in the respective dimension within which one 
half of the total energy of the i"' resolution kernel is confined under noise free conditions. 
They are calculated from the following expressions. 







r/(jr",y',z") 



Ei 



(2.29) 



//,(.r'.y',z') = 2 n^y 

V ;=i 



2 V 



r, (.v", y". r") 



£, 



(2.30) 
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(2.31) 



V(x\y\z') = 2 




rj{x".y'\ r") 



where 







(2.32) 



is the total energy of the resolution kernel and rj is the energy contained in the/' el- 
ement of the resolution kernel. In three dimensional space, the center of the box as- 
sociated with the /'* diagonal element of R is located at (x'. r', 2 '). and the center of the 
box associated with the / element of the i-'- resolution kernel is located at (.v", /',;") . 
Thus. Ax. Af. and Az are the separation distances in each dimension, respectively, be- 
tween these locations. A significant strength of measuring system performance this way 
is that the RMS error and resolution do not depend on the available data. Instead, error 
and resolution can be determined once the eigenray paths through the discretized ocean 
are determined and the covariances of the unknown variable and random error are 
known. 
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III. EIGENRAV PREDICTION AND INVERSION CODE 

DEVELOPMENT 



In this chapter, we first discuss eigenray patli prediction for the Greenland Sea array. 
Eigenrays are paths of acoustic signals extending between sources and receivers. Accu- 
rate eigenray paths are needed to establish the forward problem. Errors in ray posi- 
tioning and travel time calculations introduce errors to the inverse solution. We then 
discuss the modifications that have been made to the inversion code and their effects on 
system performance. 

A. EIGENRAV PATH PREDICTION 

The construction of the linear operator is dependent on ray path information. 
Knowing the amount of path length a ray has in each box is essential for the estimator 
to determine the distribution of unknown structures. Collectively, these structures result 
in a travel time perturbation datum as an acoustic signal passes from source to receiver. 
We use an algorithm based on ray tlieoiy to predict eigenray paths since ray theorx' 
provides a simple physical description and the equations used in modeling are straight- 
forward. 

The ray trace algorithm is a fourth order Runge-Kutta numerical integration tech- 
nique applied to Equation (2.4). The depth of each ray is calculated at 1000 range steps. 
Details of the fourth order Runge-Kutta method are given in Gerald (1989). As the ray 
reaches significant points where the vertical sense of motion reverses, i.e. turning points 
and surface reflections, the sign of the second term of Equation (3.1) reverses. (Since 
we have not included bottom intereacting eigenrays in our forward problem, we will not 
consider bottom reflections in our discussion here.) The Runge-Kutta method cannot 
be applied through a step containing a significant point. Therefore, in the vicinity of a 
turning point or a surface reflection, we apply a method which assumes a locally con- 
stant sound speed gradient at depths near the significant point. Under constant gradient 
conditions, the ray path is circular. From Ugincius( 1970), we have an expression for 
curvature of the ray path, 

K- = -y [gjr - cos 0) -l-g,,(2 cos Ov - 1)] (3.1) 
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wliere is the sound speed gradient and is the gradient of the medium speed. Tlie 
radius of curvature. X?. is just the reciprocal of k . Since we use a motionless reference 
ocean in the forward problem, the expression for the radius of curvature simplifies to 



^? = - 



c 

g^ cos 0 



(3.2) 



which shows is inversely related to the sound speed gradient with smaller radii re- 
sulting from stronger gradients. In turning point cases, the ray is projected through the 
step containing the point at a constant 5?, then the Runge-Kutta integration continues. 
In surface reflection cases, the ray is projected to the surface and continues downward 
at the same however, for the downward path, the center of the circle for 5? is sym- 
metrically relocated about a vertical line extended through the point of reflection. As in 
the turning point case, the Runge-Kutta integration continues at the completion of the 
step until the next significant point is reached. Figure 2 shows examples of ray tracing 
through steps containing significant points. 

To establish the forward problem, a reference ocean must be chosen. The associated 
sound speed profile needs to closely approximate the true conditions expected in the 
Greenland Sea so that model predicted rays can be associated with the data from the 
tomographic array, \\dien the transceiver moorings were deployed, sound speed profiles 
were generated at each site from CTD data. The profiles indicate a thin, warm layer of 
relatively fresh water lies above a sharp density gradient. The sharp gradient extends 
20-50 m below the top layer. Below the gradient, the water column approaches adiabatic 
conditions at 500-1000 m. and is essentially adiabatic below 1000 m. For simplicity in 
modeling, we treat the sound speed profile as a function of depth only. The reference 
profile is based on a profile taken near the center of the array. The profile was taken 
at the beginning of the data collection period on yearday 264 near .Mooring 6 and is 
shown in Figure 3. Yearday 1 is defined as 01 Jan 88 throughout this work. The pres- 
ence of the sharp and rapidly changing gradients in this profile near the surface in con- 
junction with the source and receiver locations proves to be challenging for our ray 
tracing algorithm. 

Under certain circumstances, small difTerences in sound speed profiles can result in 
important output difTerences from our ray tracing algorithm. The top 40 meters of two 
similar profiles with slightly dilTerent layer depths are shown in Figure 4. From 40 me- 
ters to the bottom at 3000 meters, the profiles are identical. Either profile is represen- 
tative of the conditions near Mooring 6. We have depicted the deeper layer (DL) case 
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I 

Figure 2. Path Segment Gcoiiicti y Througli Significant Points: a) 1 urning point 

example, b) Surface rcllcction example. 
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Figure 3. CTD Data at Mooring 6 on Ycarday 264: From Worchester and Howe 

(19S9). ^'earclay I = 01,01/88. 
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as a possible result of entrainment of five meters of the sharp gradient found in the 
shallower layer (SL) case by mechanical means. 

Typically, eigenrays occur in bundles of four distinct rays due to the near surface 
location of the tranceiver elements and near half channel conditions of the Greenland 
Sea. Two rays in the bundle are launched downward within a few tenths of a degree of 
each other. The deeper ray arrives at the receiver from below while the shallower ray 
reflects ofl'the surface just prior to reaching the receiver, then arrives from above. The 
other two rays in the bundle are launched upward at nearly the same magnitude launch 
angle as the downward rays. Again, the ray launched at the steeper angle arrives from 
below while the shallower ray reflects just before the receiver and arrives from above. 
Since the magnitudes of the receiving angles are nearly the same for all four rays and the 
bundles are several degrees apart, the bundles are identified by a nominal receiving an- 
gles. To identify ray paths as eigenrays, our algorithm uses a tolerance of ± 100 meters 
vertical displacement from the receiver as a threshold. This threshold usually allows us 
to identify all four eigenrays at some nominal receiving angle using 0.1° separation in 
launch angle. 

Eigenray paths between Mooring 4 and Mooring 5 with launch angles ±0 are 
shown in Figure 5 where 1 ° <0 < 14°. Figure 5a shows paths based on the DL profile. 
The SL case is shown in Figure 5b. With only a five meter difference in layer depth, 
several of the eigenrays predicted for the DL case are not predicted for the SL case. The 
rays of the SL case reach the top layer at lower angles and continue in the layer for 
longer ranges before returning to the lower ocean. The rays of the DL case reach the 
surface layer at higher angles and have less path in the high speed water, which results 
in a significantly dilTerent path from the SL case. Rays leaving a source at the same 
launch angle. 6 = —9.03°, in the two different environments are are shown in Figure 5c. 
In the DL case, the ray is projected to pass less than 100 m from the receiver and is 
considered to be an eigenray. In the SL case, a ray with the same launch angle is 
projected below the receiver by over 200 meters. 

The selection of predicted eigenray paths used for the inversion of data from the 
Greenland Sea depends on which paths best match the signal processing results of the 
array data. The data we have received for paths between Mooring 4 and .Mooring 5 
indicate, there are four distinct paths at a nominal receiving angle of ± 13° , four paths 
at ± 10°, and the slowest arrival which occurs near 0°. Therefore we have chosen to 
include the eight steepest eigenrays predicted in the DL case for our inversion of 
Greenland Sea data. These rays are at nominal angles that best match the true ocean 
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SL (solid) DL (dash) 




Figure 4. Near Surface Sound Speed Profiles: Both profiles arc representative of 

the profile at Mooring 6. 

data. In Chapter IV, we will continue our discussion of the inversion of data from the 
deployed array. 

The presence of strong and rapidly changing gradients in our reference sound speed 
profile causes another problem for our ray tracing algorithm. At relatively shallow 
launch angles, typically with 0 between + 5", the algorithm fails while projecting the ray 
through a turning point in the strong, near surface gradient. The exact cause of failure 
needs to be fully investigated, but a preliminary look indicates a higher degree of nu- 
merical accuracy is required and ray behavior in the vicinity of rajiidly changing gradi- 
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Figure 5. Sound Speed Effects on Eigenraj’ Patlis: Launcli angles: 7 ° ^ 0 ^ 14° 

a) DL case cigenrays. h) SL ease eigenrays. c) Comparison of rays 
launched at 0 = — 9.03° using DL and SL profiles. 



21 



enls is not accurately depicted. Based on these results, we can expect the accuracy of 
the ray tracing output to degrade at shallower launch angles, in the presence of high 
sound speed gradients, and in the presence of rapidly changing gradients. 

In conducting our studies involving the entire array, we select 101 eigenrays along 
the 15 paths between moorings shown in Figure 6. The geometry of the array results in 
three general types of ray patterns: 1) Short paths between outer transceivers, 2) long 
paths between outer tranceivers. and 3) paths between the outer tranceivers and the 
center transceiver. An example of each is provided in Figure 7. 

.More eigenrays were predicted, but not used for three reasons. First, we have less 
confidence in the accuracy of ray traces of shallow eigenrays. Also, initial results indi- 
cate that distinct travel time arrivals for shallower ray paths may not be available, and 
finally, we feel using these 101 eigenrays in our 500 box discretized ocean model is suf- 
ficient to determine sensitivity effects without taxing the limits of the microcomputer 
used to conduct the experiments. 

B. VERTICAL LAVERS AND RMS VARIABILITY DISTRIBUTION 

The estimator arrives at a more realistic solution if a priori information guides it 
towards a statistically sensible solution. Since we generally expect to find more strongly 
defined signatures of the unknown structures in the upper portion of the ocean volume, 
we wish to focus the resolving capability of the estimator in this region. To accomplish 
this, we make two adjustments to the inversion code. We vary the thickness of the ver- 
tical layers rather than using equally spaced layers and we describe the vertical variability 
distribution of unknown variable to the estimator through the variable's covariance 
matrix. Varying the thickness of layers allows us to use more boxes of smaller size to 
describe the region of greater variability without increasing the number of discrete boxes. 
Thus, the computational effort to generate a solution is not significantly affected. 

To show the effects of our adjustments on model performance, we use an ocean 
volume which is 240 km x 240 km in the horizontal, and 3 km deep. The horizontal di- 
mensions are evenly divided into ten 24 km segments. There are five vertical layers. 
Starting at the surface, the layers are 100 m. 200 m. 450 m, 750 m. and 1500 m thick, 
respectively. Because of the uneven vertical spacing, the lower boxes generally contain 
a larger portion of eigenray paths than if the model used evenly spaced layers. As a re- 
sult. the estimator places greater weight on the lower boxes when generating solutions. 
Figure 8 depicts the resolution kernel peak values along a vertical slice using variable 
thickness layers (VTL). Figure 9 shows the resolution kernel peak values for an equally 
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Figure 6. Plan \'icw of Array Geometry 

spaced (EQS) discretized ocean having layers of 600 m each. Both eases use the same 
eigenray patlis and the slices are taken parallel to the y-axis along the seventh increment 
in X, that is for I44km < x < 16Skm (sec Pigurc 6). 1 hese slices essentially parallel the 
ray paths between Mooring 1 and .Mooring 3. Both estimators assume a constant R.MS 
sound speed perturbation equal to 2 nrs throughout the ocean. 

The resolution kernel peak values arc a measure of the resolving capability of the 
estimator at a particular location. 1 hey reflect the way the estimator views the data 
collection process as an acoustic signal travels along an eigenray path. In the EQS ap- 
proach, the estimator assumes that a greater portion of the signal was accumulated 
above 1000 m since most of the eigenray paths are found in this region, however, this 
approach uses fewer boxes to describe the structure of the upper ocean. I he V'fL ap- 
proach, iit contrast, uses more boxes to describe the upper region and thus gives more 
detail on the solution there. However, the V I'L approach also shifts resolution towards 
the lower ocean where the system has long paths passing through the thick lowest layer 
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Figure 7. Eigciiray Pallis along Vertical Slices: a) Paths between Moorings I and 

2a. h) Patlis between Moorings I and 3. c) Paths between Moorings I 
and 6. 
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Figure 8. Resolution Kernel Peak Values (^TL approach): A value of 1.0 implies 

a bo.x is well resolved. 

of boxes. The increased detail on system pcrfonnancc in the upper ocean and the in- 
creased number of estimated points provided by using smaller boxes arc characteristics 
of the VTL approach uiiich we wish to exploit since we expect greater variability in this 
region. At the same time, we wish to limit the clicct of shifting system resolution to- 
wards the lower ocean. 

In the oceans described above, the estimators do not expect greater variability in the 
unknown field in any particular region. 1 his allows the estimator the freedom to gen- 
erate solutions which show features of unusually high intensity in regions where we re- 
alistically do not expect them. This occurs if the solution meets the requirements of the 
estimator, that is, it fits the data and provides the least mean square error. We can in- 
form the estimator to anticipate greater RMS values of the unknown field in the upper 
ocean through the covariance of the estimate. In this way the estimator expects little 
travel time perturbation will accumulate in the acoustic wave as it passes through the 
deeper ocean layers and greater variability is encountered in the upper regions. As a 
result, the resolution kernel peak values show an upward shift in model resolution, off- 
setting the downward shift caused by varying the thickness of layers. In Figure 10, we 
show a vertical slice of the results of using the VTL approach with a vertical dependent 
variability function. The vertical slice is the same as that used in f igures S and 9. In 
our example, the RMS sound speed perturbation is modeled as an exponentially decay- 
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Figure 9. Resolulioii Kernel Peak Values (EQS approacli): A value of 1.0 implies 
a box is well resolved. 

ing function with depth. An average \alue is tlien calculated for each layer. For the 
example shown in Figure 10, we have characterised the function by a surface RMS value 
of 6 in.'s and an c-folding depth of 1000 in. The average value for the entire ocean using 
this function is 1.9 m s, which closely approximates the 2 m,'s RMS Sc u.sed for the ex- 
ample shown in Figure 9. 

Though we show an improvement in resolution in the upper region of the ocean, it 
does not come without cost. As we improve resolution in a region, we also allow more 
random noise to pass through the estimator. The effect can be seen in the standard de- 
viation of the estimate, and thus a less reliable solution. The balance of higher resol- 
ution and lower random error is a choice which must be decided on the basis of 
resolvable eigenrays paths, array geometry, and the amount of experimental noise that 
can be expected. Generally, we wish to increase system resolution with a minimal in- 
crease in standard deviation of the estimate. Figure 1 1 shows the associated upward 
shift in standard deviation of the estimate corresponding to an upward shift in system 
resolution. The RMS Sc distribution used for the examples shown in Figures 11a and 
lib are the same as that used for examples shown in Figures 9 and 10, respectively. 
Typically, bias tends to dominate the total error of the system and the increase in ran- 
dom error is less important. 
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Figure 10. Resolution Kernel Pe.nk \’alues Using Depth Dependent 
^^^rial)ilily: A value of 1.0 implies a box is well resolved. 



In our modeling of the Grecnlaml Sea. we as«:ume the variability of sound .speed 
perturbations and currents decrease exponentially with depth, though each variable may 
be parameteri7ed by a dincrent folding depth. Rettconable surface R.MS dc values may 
be estimated from sea surface temperatuie measurements derived from satellite data and 
ship observations. Surface current variability may be estimated from drifting buoy and 
CTD stations at the array moorings. Since surface values may be more reliably deter- 
mined. the focus of our sensitivity study is placed on the effect of uncertainty in sub- 
surface variability. 1 he computer simulated ocean on which inversions were performed 
had a surface R.MS <)c of 6 m/s and a folding depth equal to lOOO m. The estimator is 
provided with the s;ime surface value and then we compare model performance using 
folding depths of 500 m, 1000 m. and 3000 m. A similar sensitivity study is conducted 
in current tomography. Ihe surface RMS curient for the computer generated ocean and 
the estimator is 10 cm s. The ocean wac generated using a folding depth of 700 m. 
System performance is evaluated using folding depths 350 m. 700 m, and 2100 m in the 
estimator. Tigiire 12 shows the depth dependent variability curves at different folding 
depths. ’1 he estimator and the simulateil ocean both use an average value based on 
these curves for each layer of the discretized ocean. 

A longer folding depth used in the estimator implies the system expects greater var- 
iability in the unknown field at depth. We define the half-depth, '/oZ.. as the percent 
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Figure 11. S(aiid.nid Deviation of tlie Estimate (ni/s): a) RMS cic is constant 

throughout ocean volume, b) RMS Sc is exponentially decreasing with 
depth. 
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Figure 12. Depth Dependent Vai i.Thility Curvc.s at Various Folding Deptlis: a) 

Density tomograpliy curves, b) Current tomograpliy curves. 
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of the ocean measured from the surface which contains 50% of the integrated variability 
over the ocean depth. It is expressed as 



where /is the folding depth and = —3000 m. is the bottom depth. Note that we have 
defined /so that it is positive in our application here, that is/> 0 implies variability with 
depth decreases. However. is negative to be consistent with our upward pointing 
positive 2-axis with z = 0 m at the surface. It is interesting to note that the half-depth 
depends only on the folding depth and ocean depth, but not on the surface value. This 
allows us to make comparisons between the two areas of tomography. 

In Figure 13a. we show the top three layers of a simulated sound speed perturbation 
field using the VI L approach. Figures 13b, 13c, and 13d show the corresponding esti- 
mated fields based on folding depths of 500 m, 1000 m, and 3000 m, respectively. The 
simulated field has ® oZ< = 21.5%, or half of the integrated variability is contained in 
approximately the top fifth of the ocean. An estimator using a folding depth at 500 m 
expects half of the integrated variability to be contained in roughly the top lO^'o of the 
ocean. As the acoustic signal passes through lower regions of the ocean, the signal en- 
counters higher variability than the estimator expects. The estimator interprets this 
"energetic" data as the result of more intense structures in the upper ocean than are truly 
present. The solution is constructed accordingly. The opposite effect occurs when the 
assumed folding depth is larger than the true value. With / = 3000 m. the estimator 
assumes the ocean contains half the variability in the top 40% of the ocean and tends 
to smooth the variability over greater depths than found in the true solution. 

We have made the assumption that variability is an exponential distribution. Under 
these conditions, the estimator is more sensitive to underestimating the folding depth 
than overestimating it. Figure 14 shows the relationship between /and %Z_,. The slope 
of the curve indicates system sensitivity. For a given folding depth, the slope is greater 
when the assumed value is less than the true folding depth, indicating the system's in- 
creased sensitivity to underestimation. Since the sensitivity measured by %Zj is de- 
pendent only on folding depth and ocean depth, we can also see from Figure 14 that 
when true folding depths are small, the system is more sensitive to both underestimation 
and overestimation than when the true folding depth is large. 

For current tomography, we find a similar trend. Figure 15a shows a simulated 
current field in the top three layers of the ocean. Figures 15b, 15c, and 15d show the 
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Figure 13. Peituihafion Estimator Response Using Various Folding Deptlis: a) 

Simulated Sc (m/s), / = 1000 m. b) Estimated Sc (m's),/ = 500 m. 
c) Estimated Sc (m's), / = 1000 m. d) Estimated Sc (m/s), / = 3000 
m. 
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Figure 14. Half-dep(li as a Function of Folding Depfli: Slope indicates model 

sensitivity to estimating true folding depth. 

corresponding estimates based on folding depths of 350 m, 700 m, and 2100 m, respec- 
tively. As in our density tomography, the estimator is more sensitive to underestimating 
the folding depth than overestimating. 

'I hese results of these sensitivity studies are based on the assumption that the true 
variability distribution is an e.xponentially decreasing function with depth. If the dis- 
tribution is significantly dilTerent from this, our results may not apply, however, the 
variation of the half-depth as a function of some other applicable parameter may still 
be a useful measure of .sensitivity. 
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Figure 15. Curienl Estimator Response Using Various Folding Depths: a) Simu 

lated currents, / = 700 in. b) Estimated currents, / = 350 m. c) Esti 
mated currents,/ = 700 m. d) Estimated currents,/ = 2100 m. 
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IV. PRELIMINARY INVERSIONS USING GREENLAND SEA DATA 



The discussions from the previous chapters have given us insight into the expected 
system performance in the Greenland Sea environment. However, the use of simulated 
oceans cannot explore ever}’ aspect of system performance. Only using actual data from 
the tomographic array can disclose the extent of success the system will have as an ocean 
monitoring tool. 

At present, a set of travel time measurements from transmissions between Mooring 
4 and Mooring 5 has been processed for our analysis by WHOl. Acoustic signals were 
transmitted six times a day on ever}’ third day starting in late September. 1988, and 
ending in late July, 1989. covering the winter cooling season. The complete set of data 
for transmissions from Mooring 4 to Mooring 5 is used in our preliminar}’ study. The 
quality control "goodness estimates" of some data are. however, less than optimal. In 
future studies, a criterion for the minimum acceptable goodness estimate may be set. 
For our preliminar}' study, we use for each ray an average of all six transmissions to 
determine the travel time of the day. The data indicate three bundles of eigenrays were 
resolved from signal processing at WHOl. We will make reference to each bundle by the 
magnitude of its nontinal receiving angle. The bundles consist of four rays at 1 3°, four 
rays at lO". and one ray at 0°. Using a sound speed profile based on CTD taken near 
the center mooring, our ray tracing algorithm predicts four eigenray paths near I3°and 
four paths near 9° which we use for our data inversion. Since the algorithm fails at 
shallow angles, the 0° ray path is constructed by an alternate method. We trace an a.xial 
path along the minimum sound speed axis, which is located approximately 30 m above 
the source depth. The sound speeds at the source and receiver depths are less than .3 
m,s greater than the ntinimum sound speed of the reference profile. We feel this is a 
good approximation for the 0® ray since WHOl identified this arrival by using the last 
signal peak which was ver}’ much bigger than noise (Pawlowicz.1991). To display our 
results, we have taken a vertical slice in y-range as shown in figure 1 6a. This slice con- 
tains .Mooring 4 and Mooring 5. The selected eigenray paths, projected onto the vertical 
slice, are shown in in figure 16b. 

To evaluate system performance and analyse the data, we divide the ocean volume 
into 5C)0 boxes using the VTL layering scheme described in Chapter HI. Using this 
scheme, the source at .Mooring 4 is in the top ocean layer where 48km < x < 72km, while 



34 



d«otn \mi 



• « Mooring 4 




b) 




Figure 16. Verlical Slice and Eigenray Paths for Moorings 4 and 5: a) Orientation 

of vertical slice used to display model output and system performance 
in the vicinity of Moorings 4 and 5. b) Selected eigenrays projected 
onto the vertical slice. 
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tlie receiver at Mooring 5 is in the second layer from the top where 
16Skm<-v< 192km We set the a priori knowledge for the estimator as shown in Table 



Table 3. SYSTEM PARAMETERS 



Parameter 


Assumed 

Value 


Surface RMS cic 


6 m s 


Folding Depth 


1000 m 


Random Noise 


1 ms 


Correlation Length (x.y) 


30 km 


Correlation Length tz) 


300 m 



Once the system parameters are established, the estimator's performance may be evalu- 
ated. 

Figure 17 displays the resolution kernel peak values. The resolution of the system 
is highest in the vicinity of the moorings since this is where the three sets of rays con- 
verge. The higlier resolution in the upper 1000 meters is due to the e.xponential distrib- 
ution of variability and path length contributions from both the 10° and 13° rays in this 
region. 

Another measure of system resolution is the minimum resolution length. As shown 
in figure IS. the minimum resolution length in the y direction is about 60 km near the 
top center portion of the vertical slice and near the source and receiver sites, and in- 
creases rapidly where there are few or no eigenray paths. The relatively low values are 
due to the convergence of all three ray groups in these regions. Ocean features smaller 
than the minimum resolution length are not adequately resolved and will appear 
smoothed in the estimate. That is. we expect the estimator to generate less intense and 
more spread out features than the true structure if the true solution has a length scale 
less than the resolution length at its location. Since we are using a two point system 
which has a north-south (y-direction) orientation, there is virtually no resolution of fea- 
tures in the x-direction. As more data become available from other moorings, the esti- 
mator will gain more information in all dimensions. Figure 19 shows the vertical 
minimum resolution length which also tends to be lowest where the three groups of rays 
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Figure 17. Resolution Kernel Peak Values for Greenland Sea Data Inversion: A 

value of 1.0 implies a box is well resolved. 

converge. However, it is interesting to note the unusually large minimum resolution 
length of 1700 m in the top center of the slice. It appears to be an anomalous feature 
caused by the arbitrary choice of the the coordinate system. The exact cause of this lo- 
cally high value, however, needs further investigation. 
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Figure 18. Y-Diredion Minimum Resolution Length (km) 




Figure 19. Z-Directioii Miniimini Resolulion Leiiglh (ni) 



The fact that the system has low resolving power is not surprising since there are 
only three bundles of eigenrays involved. Hy using depth dependent variability in the 
estimator, wc have concentrated the resolution of the system in the upper ocean. Since 
RMS 8c varies with depth, wc normalise the RMS error to show that the system mini- 
mizes error in the upper ocean, the region of interest. Figure 20 depicts the RMS error 
as a percent of the assumed variability in the layer. 

Producing sound speed perturbation estimates with our estimator requires wc con- 
vert the travel time measurements from the array to travel time perturbations with re- 
spect to a reference ocean. Once the sound speed perturbations have been estimated, 
we can use Equation (2.6) to determine an estimated sound speed field of the Greenland 
Sea. 

For each of the nine eigenrays identified by WIIOI, a time averaged value was cal- 
culated. The corresponding travel time difl'ercnccs between the nine time averaged 
Greenland Sea values and the nine predicted cigenray arrivals based on ray theory range 
from 200-300 milliseconds. The differences are largely due to structural difTercnccs be- 
tween the model ocean used to predict our ray paths and the "time averaged" Greenland 
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Figure 20. Percent RMS Error 

Sea, however, model errors in mooring positioning and in travel time calculations from 
our ray theory based algorithm may also contribute significantly to the difierences. 
Based on receiving angles and sound speed data taken at the mooring sites, we have a 
high degree of confidence in the positioning of the ray paths within the discretized boxes. 
Even in a worst case scenario, the 200-300 millisecond differences would imply less than 
500 meter differences in corresponding path lengths over total path lengths in excess of 
120 km. These dificrcnccs will not alfcet the linear operator of the estimator greatly. 

Errors in travel time calculations even as low as tens of milliseconds, however, may 
contribute significant errors in the estimates since they arc treated by the estimator as 
part of the data. The scale of ST in the Greenland Sea is on the order of tens of milli- 
seconds, based on differences between daily travel times and the averaged travel time of 
each of the rays. Though we cannot eliminate these errors, we can separate the errors 
from our time scries of estimates by basing the travel time perturbations on the time 
averaged ocean. By choosing the time averaged ocean as our reference state, we can 
observe the response of the estimator over the winter cooling season without superim- 
posing the effects of model travel time errors. T he structure of the time averaged ocean 
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is not known, however, and may have range dependent features of interest. We can find 
this structure with respect to the known sound speed structure used to predict our ray 
paths by inverting the travel time differences between corresponding rays of the model 
ocean and the time averaged data from the real ocean. From this information, we can 
find the sound speed structure of the Greenland Sea. 

Figure 21 shows the estimated Sc field of the time averaged ocean referenced to the 
ocean used to predict our ray paths. To produce this field, we increased the horizontal 
correlation lengths of the estimator to 100 km based on the expectation that the time 
averaged ocean is more highly correlated than any of the daily perturbation fields we 
estimate. The solution shows a reasonable trend. The time averaged ocean is expected 
to be cooler in the upper ocean since the ocean used for ray prediction is based on a 
profile taken in September, a relatively warm time in our data collection period. How- 
ever, the accuracy of the predicted structure is alTected by travel time errors from the ray 
tracing algorithm and mooring positioning in the model. An independent source of data 
is needed to verify the extent of the elTects, but positioning errors may be reduced as 
more data sets of the array are included. 

Once the system produces the first estimated perturbation field, there is little in- 
crease in computational effort to generate a time series of solutions using the same esti- 
mator. To demonstrate the response of the estimator using the Greenland Sea data from 
WHOI. we plot a time series of the average dc based on the top two model layers at the 
position of Mooring A and Mooring 5. Ihc time series are shown in Figures 22a and 
23a respectively. The general "spikincss" of the plots is the result of having a two point 
system which bases solutions on only three groups of eigenrays. It is dilTicult for the 
estimator to precisely position features with so little information. The features appear 
to move more abruptly in sequential solutions than could be physically be expected, in- 
dicating instability of the solutions. The inclusion of data from other moorings will in- 
crease the estimator's knowledge of the structures in the vicinity of Moorings 4 and 5, 
and the solutions arc expected to show greater stability over time. However, even with 
instability aflecting the short term solution, there are indications of trends at longer time 
scales detected. The seasonal cooling and warming through the fall, winter, spring and 
early summer time frame arc evident at both moorings. Since sound speed is directly 
related to temperature (Mackenzie. 1981 ), Sc decreases with winter cooling and increases 
with spring and summer warming. Events at the 20-90 day scale are extracted from the 
time series by apphing a time domain filter with cutoff frequency .05 cycles day at each 



41 




Figure 21. Sound Speed Perturbations of the Time Averaged Ocean (in/s): The 

sound speed structure of the time averaged ocean is the sum of this field 
and the sound speed structure of tlie ocean used to predicted eigenray 
paths. 

of tlie point estimates. The time series of the filtered solutions near Moorings 4 and 5 
are shown below their respective unfiltered versions in Figures 22b and 23b. 

Environmental temperature data taken at Moorings 4 and 5 sites at a depth of 100 
meters has been provided by WIlOl and are displayed in Figures 22c and 23c. Com- 
paring the filtered estimates to the respective temperature data provides evidence that 
the two point tomographic system is detecting synoptic scale events. Mooring 5 partic- 
ularly shows good correlation between temperature and the local Sc estimates. Of 'spe- 
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Figure 22. Sound Speed Estimates and Temperature Data Near Mooring 4 : a) 

Undltcred time series of estimates, b) Filtered time series of estimates 
using cutofT frequency = .05 cycles/day. c) d cinpcrature data provided 
by WIIOI. 
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Figure 23. Sound Speed Es(ima(es and Temperature Data Near Mooring 5: a) 

Unfiltcred time series of estimates, b) Filtered time series ol estimates 
using cutoff frequency = .05 cycles 'day. c) Temperature data provided 
by WIIOI. 
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cific interest are the detection at Mooring 5 of the strong cooling event near yearday 360 
and the corresponding responses to wanning and cooling from yearday 390 to about 
yearday 560. 

As an example of system output, a sequence of seven <5c fields covering a three week 
period from yeardays 462 to 480 at three day intervals is displayed in Figures 24-30. The 
temperature record at Mooring 5 indicated that a strong warming event followed by a 
strong cooling event occurred during this period while a weaker pair of warming and 
cooling events occurred near Mooring 4. The output of the estimator is consistent with 
the temperature data and reflects these trends. 




Figure 24. Sound Speed Perturbation Estimates (m/s), Yearday 462 
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Figure 25. Sound Speed Perlui halion Estimates (m/s), Yearday 465 




Figure 26. Sound Speed Perturbation Estimates (m/s), Yearday 46H 
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Figure 27. Sound Speed Perturbation Estimates (m/s), Yearday 471 




Figure 28. Sound Speed Perturbation Estimates (m/s), Yearday 474 
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Figure 29. Sound Speed Perfui b.ntion Es(iinn(es (in/s), Yeardaj' 477 




Figure 30. Sound Speed Peilurbation Estimates (m/.s), Yearday 480 
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V. CONCLUSIONS AND RECOMMENDATIONS 



In this thesis, we discussed the advantages and problems of using a ray theory based 
algorithm to establish the forward problem, developed an inversion code which could 
be applied to Greenland Sea data, and tested this code using available Greenland Sea 
data. Here we ofl'er some conclusions and recommendations for further study in this 
area. 

A. CONCLUSIONS 

First, we looked at the advantages and problems associated with using an algorithm 
based on ray theory to model eigenray paths for the Greenland Sea tomography array. 
Ray theory provides a simple depiction of eigenrays paths and the algorithm is coded 
with uncomplicated mathematics. However, CTD data from this area shows an acous- 
tically dilFrcult enviroment may exist at least during part of the year and degrade the 
reliability of the predicted eigenray paths. Strong gradients and rapidly changing gradi- 
ents in the vicinity of the transceiver elements can strongly affect model output. Rays 
which are predicted to be turned in the gradient or near the sharp change in gradient are 
strongly influenced by small changes in the sound speed profile in this region. These 
significant differences reduce the reliability of the resulting ray path. The algorithm of- 
ten fails when predicting paths for rays launched at shallow angles. Rays launched at 
steeper angles pass through these regions at greater angles and are less affected. These 
ray paths are more reliable. 

Data prepared by WHOl indicates resolvable eigenrays between Mooring 4 and 
Mooring 5 are at large enough launch angles that they are not turned in the sharp gra- 
dient. Thus, our algorithm is adequate to model the forward problem in this environ- 
ment. 

In our inversion code development, we investigated the effects of using variable 
thickness layers vice an equally spaced scheme to discretize the Greenland Sea for the 
estimator. Using thinner layers in the upper ocean and thicker layers in the lower ocean, 
we kept the total number of boxes the same. The advantages of using this scheme is that 
it provides more solution detail in the region of interest for the same computational ef- 
fort. However, since the thick lower layers contain greater path length in this layering 
scheme, system resolution shifts toward the lower ocean. However, by specifying a 
depth dependent variability in the covariance of the unknown field, we tell the estimator 
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to expect greater variability in the upper ocean. As a result of this modification, system 
resolution shifts back towards the upper ocean. For our purposes, we have modeled the 
depth dependent variability as an exponentially deceasing function with depth. If more 
accurate information is known, it could just as easily be used. 

Since the true variability may not be known, we investigated model sensitivity to 
uncertainty in variability distribution. Using an exponentially shaped variability dis- 
tribution as our basis, we looked at the eflects of underestimating and overestimating the 
characteristic folding depth of the distribution. Defining the half-depth as the percent 
of the ocean measured from the surface which contains 50% of the integrated variability, 
we have an objectiv e means of measuring sensitivity. The half-depth gives us an idea 
of how the estimator is distributing variability in the ocean. We have shown the change 
in half-depth per unit change in folding depth is always greater when the true folding 
depth is underestimated. We have also shown from the same relationship that at smaller 
true folding depths, the model is more sensitive to both underestimation and over esti- 
mation. Half-depths relationships can be developed for other variability distributions 
and applied in similar ways to measure model sensitivity. 

Our inversion method was tested by using available data from the Greenland Sea 
array as input to our model. The system used nine eigenray paths between the source 
at .Mooring 4 and the receiver at Mooring 5. The results of our preliminary application 
of the model show evidence of detection of the seasonal cooling cycle and detection of 
synoptic scale events at time scales of 20-90 days. The two point system also shows 
unrealistically large fluctuations in the solutions at periods of less than 20 days. This 
indicates instablility of the solutions, however, we expect stability to increase as more 
data becomes available from other nioorings. 

Our acoustic tomography code is ready for application to Greenland Sea data. The 
estimator can provide concentrated estimates in the region of interest and allows for a 
depth dependent variability distribution of the unknown variable. With these modifica- 
tions. our model can use tomographic data to monitor changes in ocean structure. 

B. RECOMMENDATIONS 

CTD data has shown that the near surface structure of the Greenland Sea is chal- 
lenging to model acoustically. Though our ray theory based algorithm seems adequate 
for our purposes here, the combination of the near surface gradients and transceiver lo- 
cation make this algorithm less than ideal. However, the convenience of using a ray 
theory based algorithm to model the lower ocean stills keeps it an attractive basis for 
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sound propagation prediction. We recommend that a more accurate means of estab- 
lishing the forward problem in the Greenland Sea be ijivestigated. perhaps using a hybrid 
of ray theory and normal modes. 

Based on the promising results using just a two point system, continued analysis of 
the Greenland Sea using this code and data from the remaining 14 mooring paths of the 
deployed array should yield useful estimates of ocean circulation and be able to detect 
large scale convective events. 
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